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ABSTRACT 


The  combined  mechanical  stresses  on  aircrewmen  have 
become  increasingly  acute  as  technological  developments 
extend  the  flight  envelopes  of  our  high  performance 
aircraft.  Limitations  on  the  design  of  this  type  of 
aircraft  are  frequently  dictated  by  human  tolerance.  The 
concept  of  an  analytical  model  to  evaluate  the  biomechanical 
response  of  the  human  intervertebral  joint,  under  the 
influence  of  long  term  axial  compressive  loading,  is 
important  in  assessing  the  load  carrying  capability  of 
normal  and  diseased  vertebral  segments. 

It  has  been  experimentally  demonstrated  that  healthy 
Intervertebral  joints  are  composed  of  materials  which 
exhibit  creep  characteristics.  This  investigation  is 
significant  because  it  presents  a  study  of  the  time 
dependant  behavior  involved.  An  axisymmetric  finite  element 
model  is  employed  which  incorporates  a  linear  viscoelastic 
constitutive  relation  for  the  intervertebral  disc. 
Viscoelastic  material  constants  are  found  by  matching 
one-dimensional  experimental  data  with  the  two-dimensional 
model.  Results  are  presented  depicting  displacement 
profiles  and  stress  redistributions  occurring  as  a 
consequence  of  the  inclusion  of  these  viscoelastic 
parameters  which,  for  the  first  time,  simulate  the  actual 
human  response  to  high  compressive  loads  over  a  specific 
time  span. 
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CHAPTER  i 


INTRODUCTION 


I • 1  Ba  ckg  round 

An  understanding  of  the  mechanics  of  the  intervertebral 
joint  is  of  interest  to  researchers  in  many  areas;  ranging 
from  problems  encountered  during  pilot  ejections  to  those 
associated  with  the  selection  of  suitable  disc  replacement 
materials.  At  the  present  time,  however,  even  some  of  the 
salient  features  of  the  joint's  mechanical  behavior  are 
poorly  understood  .  This  Is  a  consequence  of  the 
difficulties  in  determining  in-vivo  mechanical  properties 
and  constructing  a  realistic  analytical  model.  (Belytschko 
et.al.  (3)  ) 

Research  into  the  material  properties  of  vertebral  bone 
has  been  performed  by  Galante  et.al. (10),  Rockoff 
et.al. (20),  and  McElhaney  (17).  The  use  of  the  finite 
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element  method  to  successfully  model  the  vertebral  body  has 


been  demonstrated  by  Ba lasul ramanian  et*al.(2)  and  Hakim 
et.al. (II)* 

Brown  et*al.(6),  Nachemson  (18),  Hirsh  (14),  and 
Rolander  (21)  carried  out  experiments  to  determine  the 
static  force  deflection  properties  of  the  intervertebral 
disc  subjected  to  axial  loading*  Belytschko  et.al.(4)  and 
Spilker  (23)  employed  an  axisymmetric  finite  element  model 
to  study  the  time  independant  behavior  of  the  intervertebral 
disc. 

Meanwhile,  Kazarian  (15)  reported  creep  characteristics 
for  intervertebral  joints  subjected  to  constant  axial  loads, 
and  Kazarian  et.al.(16)  illustrated  that  the  response  could 
be  adequately  modelled  with  a  three-parameter  viscoelastic 
solid*  Using  a  one-dimensional  classical  approach.  Burns 
(7)  derived  the  material  constants  for  these  three-parameter 
so  lids  • 


1*2  Purpose 

As  a  natural  extension  of  the  work  which  has  been  done 
to  date,  this  thesis  was  undertaken  to  include  in  the  model 
of  the  in tervertebral  joint,  not  only  the  geometric 
irregularities  and  material  inhomogeneity,  but  the 
viscoelastic  creep  response  previously  noted.  The  purpose 
in  doing  so  was  to  determine  the  applicability  of  using  the 
one-dimensional  parameters  derived  by  Kazarian  et*al»(16) 
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and  Burns  (7)  in  a  three-dimensional  model.  Secondly,  to 


observe  the  influence  of  a  creeping  disc  on  the  displacement 
profiles  and  stress  redistributions  within  the  vertebral 
body  • 


1.3  An  atomy  (4) 

The  human  vertebral  column  is  a  segmented  structure  of 
24  mobile  vertebrae  separated  by  intervertebral  discs.  Each 
vertebra  consists  of  a  body  and  a  set  of  posterior  elements. 
The  center  of  the  vertebral  body  is  composed  of  soft 
trabecular  bone,  and  is  encased  circumferentially  by  a  thin 
shell  of  relatively  hard  cortical  bone.  The  upper  and  lower 
surfaces  of  the  body,  which  are  also  thin  plates  of  cortical 
bone,  constitute  the  bony  end-plates. 

The  intervertebral  disc  is  a  rather  complex  entity 
containing  both  solid  and  fluid  material.  Exterior  to  the 
disc  proper  are  two  ligaments:  the  anterior  longitudinal 
ligament  and  the  posterior  longitudinal  ligament.  The 
posterior  ligament  is  attached  dorsally,  while  the  anterior 
ligament  is  attached  ventrally,  to  the  disc. 


The  combination  of  two  vertebral  bodies  with  their 
Intermediate  ligaments  and  disc  is  what  constitutes  the 
intervertebral  joint.  During  axial  loading  the  disc  appears 
to  be  the  primary  load-carrying  structure  between  vertebrae 
(14)*  (18), (21),  and  thus  any  ligament  restraint  is  removed 
in  order  to  focus  entirely  upon  the  disc  interaction.  (See 
Fig  (1.1)  and  (1.2)) 
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1 • 4  General  Approach  and  As  sump  t lo  ns 

Initially,  analytical  techniques  were  investigated  as 
an  approach  to  modelling  the  joint*  This  investigation 
revealed  the  joint  to  be  of  such  complexity  that  classical 
methods  could  not  be  applied*  The  finite  element  method, 
however,  was  found  to  be  capable  of  handling  the  complexity, 
and  therefore  was  pursued* 

An  existing  plane-stress,  plane-strain  finite  element 
program  written  by  Hinnerichs  (12)  was  modified  to  account 
for  material  viscoelasticity  and  inhomogeneity*  Since  the 

intervertebral  joint  was  idealized  as  a  body  of  revolution, 

•\ 

the  program  was  further  altered  to  accommodate  axisymmetric 
structures*  These  modifications  were  validated  before  the 
program  was  applied  to  the  problem  at  hand* 

In  addition  to  the  assumption  of  axisymmetry,  it  was 
also  assumed  that  the  intervertebral  joint  was  symmetric 
with  respect  to  the  centerline  of  the  disc*  It  was  further 
assumed  that  the  observed  creep  was  a  quasi-static 
phenomenon;  in  which  case,  the  inertia  of  the  structure  was 
neglected*  Kazarian  (15)  showed  viscoelastic  creep  of  the 
joint  to  be  dependant  on  the  presence  of  the  disc*  That  is, 
the  vertebral  body  alone  did  not  exhibit  creep  when 
subjected  to  a  constant  axial  load*  Therefore,  it  was 
assumed  that  the  trabecular  and  compacted  bone  of  the 
vertebral  body  were  linear  elastic  materials,  whereas  the 


disc  was  idealized  as  a  homogeneous  isotropic  linear 


viscoelastic  substance*  In  other  words ,  all  of  the 
viscoelasticity  of  the  joint  was  attributed  to  the  disc. 

Based  on  these  assumptions,  a  T 1 0-T 1 1  spinal  segment 
(test  I*D.  65  Kazarian  et*al.(16)>  was  discretized  and 
studied.  (See  Fig  1.2  for  their  relative  location  in  the 


spine ) 


CHAPTER  2 

LINEAR  VISCOELASTIC  THEORY 

The  classical  theory  of  elasticity  deals  with  the 
mechanical  properties  of  elastic  solids,  for  which, 
according  to  Hooke's  law,  stress  is  always  directly 
proportional  to  strain  in  small  deformations,  but 
lndependant  of  strain  rate  Saada(22).  The  classical  theory 
of  hydrodynamics  deals  with  properties  of  viscous  fluids, 
for  which,  according  to  Newton's  law,  stress  is  always 
directly  proportional  to  the  rate  of  strain,  but  independant 
of  the  strain  itself. 

The  classical  theories  of  elasticity  and  hydrodynamics 
are  Idealizations,  and  many  materials  cannot  be  adequately 
modelled  by  one  or  the  other  of  them.  For  Instance,  a 
material  which  is  not  quite  solid  does  not  maintain  a 
constant  deformation  under  constant  stress,  but  slowly 
deforms  with  time.  When  such  a  material  is  constrained  at 
constant  deformation,  the  stress  required  to  hold  it 
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gradually  diminishes  or  relaxes*  A  material  which  is  not 
quite  liquid  may,  while  flowing  under  constant  stress,  store 
some  of  the  energy  input,  instead  of  dissipating  it.  It  may 
recover  all  of  its  deformation  when  the  stress  is  removed. 
Materials  which  exhibit  such  behavior  are  identified  as 
viscoelastic  and  can  be  adequately  modelled  by  combining  the 
classical  theories  of  elasticity  and  hydrodynamics. 


2.1  Classical  Viscoelastic  Theory 

An  extensive  study  of  the  viscoelastic  theory  has  been 
performed  by  a  number  of  researchers.  Principal  results  are 
summarized  in  texts  by  Bland(5),  Chr is t ease n ( 8 ) ,  and 
Flugge(9).  A  brief  discussion  of  the  theory  is  presented 
here  as  a  background  for  the  reader. 


2.1.1  One-Dimensional  Theory  ( 9 ) 

The  building  blocks  or  discrete  elements  of  linear 
viscoelastic  theory  are  the  spring  and  the  dashpot.  Various 
combinations  of  these  elements  serve  to  model  the  many 
different  types  of  viscoelastic  materials.  To  understand 
the  combinations  of  these  elements  one  must  investigate 
their  individual  characteristics . 

Consider  a  helical  spring  which  obeys  Hooke's  law.  Its 
constitutive  equation  may  be  written  as 
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cr=  Ee 


( 2  - 1 ) 


where  <r  is  the  stress,  c,  the  strain,  and  E  the  Modulus  of 
elasticity.  Similarly  for  the  dashpot 

(2-2) 

where  is  the  dashpot  constant  and  (#)  denotes  time 
differentiation. 

Perhaps  the  simplest  combination  of  these  elements  is 
the  Maxwell  fluid  (Fig  2.1)  in  which  the  elements  are  placed 
in  series. 


VSM- 

t 


Figure  (2.1) 


To  derive  the  constitutive  equation  for  this  model,  the 
equations  of  the  two  elements  are  written 


<r=  E* 


(2-3) 


and 


<r= 


<r~ 

e.  ~ 


(2-4) 


- ^  ^ 

It  is  required  that  the  total  strain  be  the  combination  of 


the  two  elemental  strains. 
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(T=  Ee  +  i je 


(2-9) 


and  or 


which  is  written  in  standard  form 

<J ^  £  (2-10) 

where  (in  this  case)  E  and  ^  • 

More  complicated  models  are  built  by  systematically 
combining  the  discrete  elements  with  Kelvin  solids  and/or 
Maxwell  fluids-  The  constitutive  equation  for  any  such 
model  has  the  form 

pt<r  b  £<r  +  •  «•  ^  4 —  (2-1 1  > 


Or 


If  the  differential  operators  are  defined 


then  eq  2-12  can  be  written  in  the  form 

fir  -  Q  e 

which  resembles  Hooke's  law. 


(2-12) 


(2-13) 


(2-14) 
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t 


I 


The  solution  of  Eq  2-14  for  (T  or  g,  as  a  function  of 
time  is  not  a  trivial  problem.  Perhaps  the  most  universally 
applied  solution  process  is  to  subject  the  equation  to  the 
Laplace  transformation  which  results  in  an  algebraic 
relation  between  the  transforms  ^  and  c  of  stress  and 
strain 


A 


rii- 


where  is  the  transform  variable. 

If  the  polynomial  operators  are  defined 


e-t 


then  Eq  2-15  becomes 


fc- 


j 


To%* 


■As 


(2-15) 


(2-16) 


(~?<r 


(2-17) 


Again,  note  the  similarity  to  Hooke's  law. 


2.1.2  Three-Dimensional  Theory  ( 9  ) 

To  this  point  in  the  presentation,  consideration  has 
been  given  only  to  uniaxial  states  of  stress  and  strain  and 
an  expression  for  the  constitutive  equation  of  a  one- 
dimensional  linear  viscoelastic  material  has  been  derived. 
Next,  a  multlaxlal  state  of  stress  is  considered  and  a 
general  constitutive  equation  for  a  three-dimensional 
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viscoelastic  material  is  found*  In  order  to  do  this,  the 
elastic  constitutive  equations  must  first  be  presented* 

In  elasticity,  the  state  of  stress  and  strain  in  a 
material  can  be  described  with  second-order  cartesian 
tensors  and  (22)*  If  the  material  is  isotropic,  then 
only  two  material  constants  are  required  to  express  the 
constitutive  equation*  These  constants  are  Young's  modulus 
and  Poisson's  ratio*  The  constitutive  equation  may  be 
expressed 


S  ■  |WW 


(2-18) 


where  [A]  is  dependant  on  Poisson's  ratio  V  • 

Another  way  of  expressing  this  relationship  is  by 
separating  the  stress  tensor  into  its  spherical  and 
devlatoric  components  and  writing 


s  =  3Ke 


(2-19) 


s  *2&c 


(2-20) 


where  K  and  G  are  the  bulk  and  shear  moduli  ,  s  and  e  the 
spherical  stress  and  strain,  s  and  e  the  deviatoric  stress 
and  strain  respectively*  If  an  elastic  material  were 
subjected  to  a  hydrostatic  state  of  stress  then  there  would 
be  a  volume  change  but  no  change  in  shape*  In  other  words, 
the  spherical  components  (Eq  2-19)  would  be  affected,  but 
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the  deviacoric  components  (Eq  2-20)  would  not. 

Similarly*  if  a  viscoelastic  material  is  isotropic*  a 
hydrostatic  stress  must  produce  a  volume  change  and  no 
distortion.  The  quantities  s  and  e  must  therefore  be 
connected  by  a  relation  such  as  eq  2-14*  and  should  be 
written 


o  r 

P  s  -  Q  e 


(2-21) 


(2-22) 


Similarly  for  the  shear  components 
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Or 


(2-23) 


(2-24) 


P  *  A 

,  Q,  r ,  Q*  which  describe  the 
viscoelastic  behavior  of  the  material  are  entirely 
independent  of  one  another. 

The  elastic  solid  is  a  limiting  case  of  the 
viscoelastic  material.  It  is  seen  that  in  this  case  the 
four  operators  are  simply  muli tplicative  constants 
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(2-25) 


p-  i  ,  Q  =  3<  ,  p'-  f  ,  Q  - 26 

In  an  elastic  solid  under  constant  load,  nothing 
depends  on  time.  In  a  viscoelastic  material,  all  stresses, 
strains,  and  displacements  which  are  material  dependant  are 
also  time  dependant.  If  the  Laplace  transformation  is 
applied  to  eqs  2-22  and  2-24  they  become 


and 


/  jv 

* 


(2-26) 


(2-27) 


These  are  algebraic  relations  which  become  identical 
with  their  elastic  counterparts  (Eqs  2-19  and  2-20)  if  the 
following  substitutions  are  made: 


This  leads  to  a  most  general  correspondence  principle: 
(5)  If  the  solution  of  an  elastic  problem  is  known,  the 
Laplace  transformation  of  the  solution  to  the  corresponding 
viscoelastic  problem  may  be  found  by  replacing  the  elastic 
constants  K  and  G  by  quotients  of  polynomial  operators  and 
the  actual  loads  by  their  Laplace  transforms. 

In  most  cases  the  solutions  of  elastic  problems  are  not 
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written  in  terms  of  K  and  G,  but  rather  in  terms  of  E  and 


V  • 


The  elastic  relations  between  E, 


r  ±}CG 

3K+& 


^  ,  K  and  G  are 


££  -  2& 

6K.  +-2G 


(22) 


(2-29) 


If  the  corresponding  viscoelastic  polynomial  operators 
for  K  and  G  are  substituted  into  Eq  2-29,  the  following 
expressions  are  obtained: 


£ 


and 

2@'£+sL<3'" 


(2-30) 


(2-31) 


This  correspondence  principle  is  of  sweeping 
generality,  but  it  does  have  its  limitations.  It  requires 
that  one  has  closed  form  solutions  to  the  elastic  problem, 
which,  for  nonhomogeneous  structures  of  complicated 
geometry,  are  unattainable. 

The  viscoelastic  creep  of  the  human  intervertebral 
Joint  is  such  a  problem.  Thus,  the  classical  viscoelastic 
correspondence  principle  must  be  abandoned  as  a  solution 
technique  and  another  method  must  be  tried,  namely  the 
finite  element  method. 
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2.2  The  Finite  Element  Method 


Numerous  authors  have  presented  works  which  detail  the 
solution  of  linear  and  nonlinear  elastic  problems  by  the 
finite  element  method.  Among  them  are  Bathe  et.al.(3), 
0den(19),  and  Z i enki ewi c z ( 2 5 ) •  Appendix  A  gives  a  brief 
overview  of  its  application  to  an  elastic,  axi-symme t r ic 


structure 

under 

initial  strain. 

The 

solution 

of 

the 

viscoelastic 

problem  by  this  technique 

is  based 

on 

the 

assump  t ion 

that 

the  material  of 

the 

s  t  rue  tur e 

can 

be 

adequately  modelled  by  a  finite  number  of  Kelvin  solids  or 
Maxwell  fluids  coupled  with  an  elastic  spring  [Adey 
et.al.(l),  White(24),  Zienkiewicz  et.al.(26)].  This 
assumption  is  satisfied  as  long  as  the  coefficient  of  the 
zeroth-order  derivative  of  eqns  2-12,  2-21,  or  2-23  is 
nonzero • 


2.2.1  One-Dimensional  Problems  (26) 

For  simplicity  of  presentation,  it  is  assumed  that  the 
structure  In  question  is  one-dimensional  and  can  be  modelled 
as  a  three-parameter  solid,  le«  One  Kelvin  solid  coupled 
with  an  elastic  spring.  The  Instantaneous  elastic 
deformation  is  dependant  only  on  the  spring  and  therefore 
can  be  separated  from  the  creep  deformation  and  solved  by 
the  finite  element  method.  The  Kelvin  solid  (Fig  2.2)  whose 


constitutive  relationship  is  given  by  eqn  2-9,  is  next 


considered 


This  equation  may  be  rewritten 


(2  —  32  ) 


which  shows  that  the  creep  strain  rate  is  a  function  of 
current  stress  level  and  total  accumulated  creep  strain* 
This  expression  is  integrated  by  Euler's  time  derivative 
method  [Ho r nbeck ( 1 3 ) ]  to  yield 


A€ 


-  _L(<r-  At. 

C  ^ 


(2-33) 


The  quantity  4tcrepresents  the  change  in  strain  due  to 
creep  at  the  end  of  a  small  time  interval  •  This  can  be 
considered  as  an  initial  strain  imposed.  The  elastic 
problem  is  solved  again  for  a  new  stress  system.  Time  is 
Incremented  and  the  whole  process  is  repeated  in  this 
incremental  fashion  until  a  maximum  time  is  reached*  The 
foregoing  algorithm  is  more  completely  discussed  in  Appendix 
A. 

In  general,  for  any  viscoelastic  material  model,  eqn 
2-14  must  be  employed,  which  can  be  rewritten 


€«.  -  q  (2-34) 

Advantage  is  taken  of  the  similarity  of  this  equation  with 
Hooke's  law.  That  is,  for  uniaxial  conditions,  the  quantity 
Va  of  a  simple  elastic  case  is  replaced  by  the  operator 
p/a  for  the  viscoelastic  case. 


2.2.2  Th  ree-Dlmenstonal  Problems  (26) 


To  generalize  this  technique  to  the  multiaxial  case  an 
appeal  is  made  to  the  multiaxial  elastic  relationships  2-19 
and  2-20  as  well  as  the  similar  viscoelastic  relationships 
2-22  and  2-24. 

From  these  relations  it  can  be  seen  that  and 

i/2G  of  the  simple  elastic  case  are  replaced  by  the 
operators  P/a  and  P/a  for  the  viscoelastic  case. 

Again  for  simplicity  of  presentation,  it  is  assumed 
that  the  material  is  elastic  in  volume  change  and  a  three 
parameter  solid  in  distortion.  With  these  assumptions,  the 
operators  become 


1/3  K-  =  P?&  =  V3K 


(2-35) 


and 

1/26  =  P '/o'  -  */(<).  *  %  °)  (2-36) 

where 

D  =  J/i  t  <2-37) 

As  In  the  one-dimensional  case  it  is  necessary  to 
separate  the  elastic  terms  from  the  viscoelastic  creep 
terms.  This  is  done  simply  by  noting  which  terms  are  time 


21 


dependant  and  which  are  not.  The  volume  change  is  assumed 
to  be  all  elastic  so  it  may  easily  be  written 

1/3K=  (  +  (0)c  (2-38) 

where  the  subscripts  "e"  and  "c"  denote  elastic  and  creep 
respectively.  The  elastic  part  is  given  by 


CV3K.)e  =  V3K.  (2-39) 

and  the  creep  part  is  given  by 

0/3f0c  =  0  (2-40) 

From  the  distortion  equation  the  parts  are  separated 
similarly.  The  elastic  part  is 


C*/2fi)e  =  ^ 

and  the  creep  part  is 


(2-41) 


(V 2&)e  =  V(%.+  %0)  (2-42) 

Now »  since  it  is  desired  that  these  relations  be 
applied  to  the  finite  element  method,  the  values  y  and 
Elastic  modulus  must  be  solved  for.  In  order  to  do  this, 
the  elastic  parts  are  grouped  together 


(tys0«  ;  V3K  ;  0/26  )e  =  ^ 


(2-43) 


Elastic  modulus  and  V  are  found  using  eqn  2-29 


*7 1C 

fiK'j.'*'  < 

and 


Cv)t 


&K\  +  1 


the  creep  parts  are  grouped  together 


(2-44) 


(2-45) 


(i/3fc.\  =  0  j  0/2G)c  =  V(p+ %!>)  (2-46) 

and  solved  for 

Cv)t  =  |/2  ;  (Ve)<.  =/6(  (2-47, 

In  order  to  use  these  quantities  in  the  finite  element 
method  eqn  2-18  must  be  returned  to*  For  the  elastic  part 
the  values  of  and  E^are  substituted  and  the  instantaneous 
elastic  deformation  is  solved  using  the  finite  element 
method*  For  the  creep  part,  the  values  for  and  E^are 
substituted  to  obtain 

=  2M %p)  (2-48) 

which  is  rewritten 

%[**}  <2-‘,» 

This  equation  Is  Integrated  by  Euler's  method  to  yield 
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[f£cl =  ~  J  <2'5(,, 

Again,  as  in  the  uniaxial  case,  the  becomes  an  initial 
strain  imposed  and  the  finite  element  method  proceeds 
similarly.  (If  the  reader  is  specifically  interested  in  the 
use  of  Eq  2-50  within  the  finite  element  formulation. 
Appendix  A  may  be  consulted.) 
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CHAPTER  3 


PROGRAM  VALIDATION 


As  was  mentioned  in  the  introduction,  numerous  changes 
were  made  to  an  existing  finite  element  program.  The  first 
step  following  these  changes  was  to  test  the  program  against 
some  known  exact  viscoelastic  solutions* 


3-1  One-Dimensional  Axial  Rod 

The  first  test  was  conducted  to  demonstrate  the 
validity  of  the  program's  application  to  a  simple  one- 
dimensional  problem*  A  homogeneous  circular  rod  having 
three-parameter  viscoelasticity  was  subjected  to  a  constant 
axial  load*  Pig*  (3*1)  shows  the  rod  and  the  finite  element 
mesh  used  to  model  it,  while  a  comparison  of  the 
displacements  with  the  exact  solution  is  presented  in  Fig. 
(3*2)*  The  depicted  solution  was  obtained  by  using  an 
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initial  timestep  of  ten  seconds  which  varied  to  a  value 
47  seconds  at  the  conclusion  of  the  problem, 
approximate  values  were  within  1%  of  the  exact  solution. 


Figure  3*1  Homogeneous  Circular  Rod 


Figure  3.2  Axial  Rod  Displacement  vs  Time 
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3 «  2  T h i ck-Wa lied  Cylinder  Subjected  t o  Internal  Pressure 


Whereas  the  first  example  tested  the  program  for 
accuracy  In  handling  a  very  simple  structure  under  an  axial 
load,  the  second  example  was  designed  to  test  the  method  on 
a  relatively  complex  problem  under  radial  loading.  The 
thick-walled  cylinder  problem  was  chosen,  because  an  exact 
solution  could  be  derived  using  the  general  correspondence 
principle  and  it  also  provided  the  desired  complexity*  Pig* 
(3.4)  shows  the  geometry  and  boundary  conditions  of  the 
problem • 


3*2*1  Exact  Solution 


In  order  to  apply  the  general  correspondence  principle, 
the  elastic  solution  was  needed.  The  plane-strain 
assumption  was  made  which  led  to  the  following  elastic 
solution  for  the  problem: 


(3-1) 


o+v)  P)L_r 

E(Al-blJ  L 


( J  -  Zv^  r* 


(3-2) 


It  was  noted  that  and  are  independant  of  the 

material  constants,  and  therefore,  in  the  viscoelastic 
problem  these  stresses  should  be  time  invariant*  The 
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however,  are 


displacement  ur  and  the  other  stress  ^  , 
material  dependant  and  therefore  should  be  time  dependant  In 
the  viscoelastic  problem*  These  quantities  were  found  by 
applying  the  correspondence  principle  to  eqs  3-1  and  3-2 


-  zpfr  Qs“-  i'e"i 


(3-3) 


U 


(3-4) 


These  relations  were  next  applied  to  specific 
materials.  It  was  postulated  that  the  volume  change  was 
purely  elastic,  while  the  distortion  behaved  like  a 
three-parameter  solid. 


VsAA 


Figure  (3.3) 


This  assumption  led  to  the  following  polynomial 
operators : 

Q'*1,  ;  <9=  (3-5) 

which  were  substituted  into  eqs  3-3  and  3-4* 

After  some  lengthy  arithmetic  and  transform  inversions 
the  following  expressions  were  found: 
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<r£  - 


2p£. 

<f-b% 


(3-7) 


where 


6K  +(j,, 


e  -  +/*, 


(3-8) 


Appendix  B  contains  a  full  derivation  of  these  equations. 


3.2*2  Finite  Element  So lut ion 

In  order  to  apply  the  finite  element  method,  a  slice 
was  taken  from  the  cylinder  and  discretized  as  shown  in  Fig 
(3.4).  The  material  constants  were  found  using  the  method 
detailed  in  section  2.2.2  and  the  problem  was  solved.  The 
finite  element  solution  was  obtained  using  an  initial 
tlmestep  of  0.01  seconds  which  increased  to  a  value  of  0.1 
seconds  at  the  conclusion  of  the  problem.  Figs  (3.5) 
through  (3.9)  present  a  comparison  of  the  finite  element 
results  with  the  exact  solutions.  In  each  instance  the 
variations  from  the  exact  were  virtually  nil.  The 
worst-case  variation  was  less  than  it 
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Exact  Solution 


Figure  3.5  Thick-walled  Cylinder  u  vs  t  r-2  in 
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Figure  3.7  Thick-walled  Cylinder  u  vs  t  r*4  In 
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Figure  3.9  Thick-walled  Cy  1  lnder  vs  t 
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CHAPTER  4 


APPLICATION  TO  THE  INTERVERTEBRAL  JO INT 


In  the  experiments  performed  by  Kazarian  et.al.  (16), 
several  intervertebral  joint  specimens  were  studied.  These 
specimens  were  made  by  cutting  the  vertebral  bodies  parallel 
to  the  lines  of  demarcation  between  the  discs  and  the  bony 
end-plates  at  levels  of  maximum  vertebral  waisting.  The 
specimens  were  then  measured  and  placed  in  the  test 
apparatus.  Each  was  subjected  to  a  dead  weight  load  of  18.0 
Kp  (Klloponds)  for  a  period  of  time  and  the  deflection  of 
the  top  edge  was  recorded. 

For  the  purpose  of  this  finite  element  study,  specimen 
number  65  (a  T10-T11  segment)  was  used.  So  that  comparisons 
with  the  experimental  data  could  be  made,  the  same  18.0  Kp 
load  was  applied  to  the  finite  element  model  as  was  applied 
to  the  test  specimen. 
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*+  •  i  Me  sh  Size  and  Material  Constant  Determination 

In  order  to  apply  the  finite  element  method  to  the 
intervertebral  joint,  it  had  to  be  discretized  in  such  a  way 
as  to  maximize  the  accuracy  of  the  solution  while  minimizing 
the  computer  resources  expended.  To  accomplish  this,  it  was 
necessary  to  systematically  reduce  the  mesh  sizes  until  the 
stresses  within  the  joint  were  observed  to  vary  smoothly 
from  one  region  to  another.  Additionally,  since  a  major 
interest  of  this  study  was  to  observe  the  bone-disc 
interface,  the  bony  end-plate  region  was  used  as  a  focal 
point  during  the  mesh  size  reduction  process.  When  the 
stresses  within  the  bony  end-plate  were  observed  to  vary  by 
15Z  from  those  of  the  previous  mesh  size,  the  model  was 
considered  to  be  refined  enough  to  capture  the  trends  of 
stress  redistribution  with  time.  Figure  4.1  shows  the 
Initial  mesh  used  in  the  study. 

Material  properties  for  the  trabecular  and  cortical 
bone  were  taken  to  be:  E-  750Kp/Sqcm  V  -  0.25  and  E- 

161 OOOKp/Sqcm  «  0.25,  respectively  (4).  Assignment  of  the 
viscoelastic  material  constants  to  the  disc  was  done  based 
on  the  properties  derived  in  Ref.  (16)  i.e.  E-  1 6 7 . 2Kp / Sq cm , 
-  64. 6Kp/Sqcra,  and  ^  -  1 1 06000Kp-Sec /Sqcm.  (See  Fig  4.2 

Ref.  (16)  fit  to  the  experimental  data)  A  problem  in 
assigning  these  constants  was  immediately  encountered  due  to 
the  assumptions  made  in  their  derivation.  Those  assumptions 
being  that  the  Joint  was  a  homogeneous  one-dimensional  rod. 
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1 


Since  the  objective  of  this  thesis  was  to  represent  the 
joint  as  a  no nhomoge neou s ,  axisymmetric  structure,  it  was 
necessary  to  deviate  from  these  constants. 

The  problem  of  inhoraoge nei ty  was  approached  first.  If 
one  has  two  elastic  rods  loaded  as  shown  in  Fig.  4.3,  the 
first  being  made  of  a  single  material,  while  the  other  made 
of  two  materials,  it  is  readily  shown  that  the  displacement 


u(D 


SJL 

AE 


for  the  homogeneous  rod  and 


(4-i  ) 


uMl 


pq-a-1  .  ik. 
A  El  At, 


(4-2) 


for  the  nonhomogeneous  rod.  The  value  of  E  as  given  in  ref. 
(16)  was  derived,  in  effect,  by  loading  the  specimen, 
observing  the  instantaneous  displacement  u(i)  ,  and  then 
applying  eqn  4-1. 


E= 


Pi 

AuU) 


(4-3) 


This  approach  was  extended  to  the  nonhomogeneous  rod  by 
equating  the  right-hand  sides  of  eqs  4-1  and  4-2. 


PJ  .  Btil  . 

A  E  "  AU  At, 


(4-4) 


Here,  all  the  values  were  known  with  the  exception  of  E^  , 


which  was  found  by  simple  algebra  to  be 


m 


gt  EE  a. 

jf(Ea-E)+*-E 


(4-5) 


Since  the  vertebral  body  was 


elastic,  the  constants 


a  and  q 

>  l 

proportioned  by  to  obtain 


assumed  to  be  linearly 
for  the  disc  were  simply 


2  3.34  Kp / Sq cm 


5 


I 


400,000  Kp-Sec/Sqcm 


Poisson's  ratio  for  the  disc  was  taken  to  be  0.48  (4). 

These  material  constants  were  assigned,  and  a  computer 
run  was  made.  The  instantaneous  displacement  of  the  top 
edge  of  the  joint  at  the  line  of  rotational  symmetry  was 
observed  to  be  0.003  cm,  which  varied  significantly  from  the 
experimental  value  of  0.014  cm.  Furthermore,  as  time 
proceeded,  the  displacement  remained  nearly  constant.  At 
the  end  of  167  minutes  simulation  time,  a  value  of  only 
0.004  cm  was  reached,  while  the  experimental  deflection  was 
0.023  cm  at  that  time.  These  large  variations  from  the 
experimentally  observed  data  were  concluded  to  be  a 
consequence  of  having  determined  E,  cj^a»  and  for  the 
disc  based  on  the  assumption  that  the  Joint  was  a 
one-dimensional  rod. 

Before  the  mesh  size  reduction  could  proceed,  the 
values  of  E,  ^p*  anc*  ^ ,  bad  to  be  found  which  would  more 
closely  match  the  experimental  results.  It  was  noted  that 
the  Instantaneous  displacement  was  dependant  on  the  value  of 
E  but  lndependant  of  c^c  and  C^y  •  The  first  step, 
therefore,  was  to  vary  the  E  value  until  the  model  matched 
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the  experimental  displacement  at  t*0.  This  was  done  quite 
readily  using  a  trial  and  error  method.  The  value  of  E 
which  best  represented  the  instantaneous  deflection,  was 
found  to  be  7.5  Kp/Sqcra.  Fig  4.4  shows  the  displacement  vs 
time  plot  based  on  this  E  value. 


The  values  of 


and  \y 


were  found  by  a  similar 


slope  of  the  displacement  vs  time  curve,  while  the  value 

had  more  influence  on  the  displacement  limit  approached  as 
time  increased.  Fig  4.5  and  4.6  show  some  intermediate 
results  obtained  by  varying  ^  and  .  The  values  for 
and  which  best  represented  the  data  were  found  to  be 

Kp/Sqcm  and  =28,  150  Kp-Sec/Sqcm. 

With  these  material  constants,  the  mesh  size  reduction 
process  was  continued .  Figures  4# 7  and  4*8  show  some 
intermediate  mesh  sizes.  Whereas,  at  each  reduction,  the 
instantaneous  deflections  varied  little  from  those  obtained 
with  the  course  mesh,  those  observed  as  time  progressed 
began  to  deviate  from  the  experimental  data.  These 

observations  revealed  that  the  value  of  E  which  was  derived 
from  the  course  mesh  model  was  good,  but  the  values  of 
and  ^  needed  more  refinement.  Therefore,  as  the  mesh  size 
wa 8  further  reduced,  minor  adjustments  to  the  ^  and  ^ % 
values  were  made.  This  process  continued  until  the  /%  and 
j  values  showed  less  model  dependance,  and  the  stresses 
within  the  joint  varied  smoothly  from  one  region  to  another. 
Furthermore,  the  stresses  in  the  bony  end-plate  were 
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observed  to  vary  by  less  than  15%  from  those  of  the  previous 


mesh  size*  Figure  4.9  presents  the  final  mesh  size,  while 
Fig  4.10  shows  the  displacement  vs  time  plot  based  on  the 
final  values  of  ^-.01  Kp/Sqcm  and  -27000  Kp-Sec/Sqcm. 


4.2  Displacement  Profiles 

In  order  to  visualize  the  displacements  throughout  the 
joint  as  time  increases,  plots  of  the  deformed  shape  were 
made  at  convenient  intervals*  Figure  4.11a  presents  the 
undeformed  shape,  while  Figs  4.11b,c,d  show  the  shapes  at 
time-0.0,  80.0,  and  160.0  minutes  respectively.  In  each  of 
the  figures,  the  deformations  were  magnified  by  a  factor  of 
2  so  that  they  could  be  more  readily  observed. 

It  was  noted  that  because  of  the  large  disparity 
between  the  elastic  moduli  of  the  vertebral  body  and  the 
d 1 8  c ,  the  greatest  Instantaneous  deformation  presented 
itself  in  the  disc.  As  time  progressed,  this  condition 
became  even  more  pronounced  in  support  of  the  assumption 
that  all  viscoelastic  effects  were  contributed  to  the  disc. 

Note  that  the  disc  appears  to  be  "squeezed"  outwardly 
in  the  radial  direction.  As  it  moves  outward  in  time,  one 
would  expect  that  there  would  be  an  accompanying  Increase  In 
shear  stress  along  the  disc-bone  interface.  Additionally, 
since  the  disc  is  constrained  from  moving  radially  at  the 
line  of  rotational  symmetry,  one  would  expect  the  radial 
component  of  direct  stress  to  show  an  Increase  in  tension  as 


the  disc  displaces  outwardly*  This  outward  movement  also 


appears  to  cause  a  "twisting**  of  the  vertebral  body.  That 
is,  the  wall  of  the  cortical  bone  on  the  exterior  of  the 
body  is  being  displaced  outwardly  at  the  disc-bone 
interface,  while  the  top  edge  is  displacing  downward  and 
inward.  One  would  expect  this  "twisting"  to  cause  the 
trabecular  region  to  be  placed  in  a  state  of  compression 
that  would  increase  with  time. 

In  the  following  section,  the  actual  stress 
redistributions  are  presented  and  appear  for  the  most  part 
to  follow  the  stresses  one  would  expect  from  intuitive 
arguments • 


4.3  Stress  Redistributions 


The  variation  of  shear  stress  through  the  joint  was  of 
interest  in  determining  the  final  mesh  size.  Vertical 
profiles  of  shear  were  taken  at  two  radii;  r-0.625  cm  and 
r«1.5  cm,  and  plotted  to  observe  the  variations.  These 
profiles  are  presented  for  t-0.0  in  Figures  4.12  and  4.13. 
It  was  noted  that  the  bony  end-plate  served  as  a  boundary  in 
which  the  shear  in  the  two  relatively  soft  adjacent  regions 
underwent  a  sign  change. 

In  addition  to  these  vertical  profiles,  horizontal 
profiles  of  all  the  stress  components  were  taken  at  three 
different  joint  levels:  within  the  disc,  the  bony  end-plate, 
and  the  trabecular  bone  region  just  superior  to  the  bony 
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end-pLate.  Figures  4.14  through  4.25  present  these  profiles 
for  times  varying  from  0.0  to  160*0  minutes.  From  these 
stress  profiles,  a  number  of  observations  were  made. 

Figure  4.14  reveals  that  the  radial  component  of  stress 
in  the  disc  increased  with  time  and  at  a  radius  of  2.2  cm, 
the  increase  after  160.0  minutes  was  approximately  30 
This  Increase  was  also  noted  in  the  circumferential  or  hoop 
stress  component  (Fig  4.16).  The  axial  stress  component 
(Fig  4.15)  was  seen  to  increase  by  about  7%  overall  during 
the  160.0  minute  interval.  The  shear  stress  (Fig  4.17) 
showed  an  interesting  trend.  For  radial  values  from  0.0  to 
about  1.75  cm,  the  shear  increased  with  time,  while  for 
greater  radial  values,  it  tended  to  decrease  or  relax.  This 
trend  to  relax  was  not  predicted  by  intuitive  arguments,  but 
is  concluded  to  be  a  result  of  the  disc-bone  interface 
geometry.  That  is,  because  the  interface  is  a  curved 
surface  rather  than  a  flat  one,  intuition  fails  to 
accurately  predict  the  stress  redist ribut ion. 

In  the  bony  end-plate,  where  the  greatest  interest  of 
this  study  was  centered,  the  first  observation  one  makes  is 
that  in  each  case,  the  stress  components  are  of  a  much 
larger  magnitude  than  for  any  other  joint  region.  This 
points  to  the  fact  that  the  disc-bone  interface  is  an  area 
which  serves  an  important  role  in  the  overall  integrity  of 
the  human  intervertebral  joint.  Furthermore,  as  time 
progressed,  the  stresses  underwent  some  significant  changes, 
which  reveal  the  region  to  be  very  active  during  creep.  The 
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radial  component  (Fig  4.18)  Increased  by  about  8%  overall 
during  the  simulation*  In  the  shear  stress  component  (Fig 
4.21)  there  was  an  overall  tendancy  to  increase  in  absolute 
value*  This  trend  was  most  pronounced  at  r~2.0  cm  where  the 
Increase  was  about  20%*  The  hoop  stress  (Fig  4.20) 
Increased  for  radial  values  from  0*0  to  1*5  cm  and  relaxed 
for  radial  values  greater  than  1*5  cm*  The  stress  reduction 
was  most  pronounced  at  a  radius  of  2.2  cm  where  a  9%  change 
was  observed.  Figure  4.19,  which  depicts  the  axial  stress 
component,  shows  the  bony  end-plate  to  be  in  a  state  of 
compression  axially.  This  component  showed  a  tendancy  to 
relax  with  time,  with  the  reduction  on  the  order  of  9% 
overal 1 • 

In  the  trabecular  bone  region,  the  stresses  (Figs 
4.22-4.25)  remained  virtually  constant.  It  should  be  noted 
that  the  final  data  point  on  these  plots  near  r-2.40  cm  is 
In  the  region  of  the  vertebral  body  which  is  composed  of 
cortical  bone.  In  that  region  (as  was  the  case  with  the 
bony  end-plate)  the  stresses  did  vary  with  time. 
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Figure  4*3  Two  Elastic  Rods 
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Figure  4.9  Final  Mesh  Size 
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Figure  4.10  Displacement  vs  Time  Final  Results 
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Figure  4.11b  Deformed  Shape  t«0+ 
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Figure  4.11c  Deformed  Shape  t-80  min< 
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Flgure  4.17  RZ  Stress/Input  P--disc 
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Figure  4*20  Theta  Stress/Input  P — Bony  End-plate 
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Figure  4.21  RZ  Stress/Input  P--Bony  End-plate 
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Figure  4.23  Z  Stress/Input  P — Trabecular  Bone 
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CHAPTER  5 


CONCLUSIONS 


Whereas  attempts  to  model  the  viscoelastic  creep  of  the 
human  intervertebral  joint  as  a  one-dimensional  viscoelastic 
structure  have  proven  successful  in  portraying  the 
externally  observed  response,  they  have  fallen  short  of 
showing  the  internal  response*  Furthermore,  when  attempts 
are  made  to  apply  the  viscoelastic  parameters  derived  from 
such  a  model  to  a  more  complex  two-dimensional  model,  even 
the  external  response  is  not  accurately  reproduced. 

This  observation  led  the  author  to  deviate  from  the  use 
of  the  one-dimensionally  derived  parameters  and  seek  to 
derive  constants  which  would  more  accurately  reproduce  the 
externally  observed  response.  In  doing  so,  it  was  found 
that  the  material  constants  based  on  the  ax i-symme t r ic  model 
of  the  Joint  varied  substantially  from  those  derived  from 
the  one-dimensional  model.  These  significant  variations  led 
to  the  conclusion  that  for  purposes  of  stress  analysis,  the 
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one-dimensional  parameters  are  Inadequate.  Stated 
alternatively,  the  intervertebral  joint  is  not  accurately 
modelled  as  a  nonhomoge neou s  viscoelastic  axial  rod. 

As  was  pointed  out  in  the  preceeding  chapter,  the 
viscoelastic  creep  of  the  human  intervertebral  disc  causes 
some  significant  stress  redistributions  throughout  the 
entire  intervertebral  joint.  These  changes  were  most 
pronounced  in  the  disc-bone  interface  region  ,  where  the 
greatest  activity  was  noted.  These  redistributions  led  to 
the  conclusion  that  time  is  an  important  variable  which 
should  be  included  when  modelling  the  human  intervertebral 
joint. 
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APPENDIX  A 


In  the  following  sections,  the  basic  concepts  and 
equations  used  in  the  finite  element  analysare  of  elastic 
materials  is  briefly  reviewed*  These  equations  provide 
background  for  the  viscoelastic  finite  element  computer 
program  developed.  The  final  section  serves  to  show  the 
incremental  nature  of  the  viscoelastic  problem  and  ties  the 
viscoelastic  constitutive  relation  into  the  elastic  solution 
method  • 

The  basic  philosophy  of  the  finite  element  method  (25) 
is  that  an  approximate  solution  to  a  complicated  problem  can 
be  obtained  by  subdividing  the  region  of  interest  into  a 
finite  number  of  discrete  elements  and  then  choosing 
appropriate  relatively  simple  functions  to  represent  the 
solution  within  each  element*  These  functions  are  simple 
compared  to  the  so-called  "exact”  solutions  which  account 
for  the  entire  region  of  interest.  In  this  section  the 
equations  associated  with  representing  an  axisymmetric  body 
as  a  finite  number  of  elements  are  presented.  The 
displacements  in  each  element  are  expressed  as  a  simple 
polynomial  and  the  equations  relating  displacements  to 
applied  loading  are  given* 

DISPLACEMENT  MODEL 

The  displacement  function  used  in  the  displacement 
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formulation  is  generally  selected  as  a  polynomial.  The 
polynomial  expression  allows  for  simple  differentiation  and 
integration.  Also,  as  the  element  size  becomes  small,  the 
polynomial  expression  permits  a  simple  approximation  to  the 
exact  solution.  A  polynomial  of  infinite  order  corresponds 
to  an  exact  solution.  However,  for  practical  purposes  the 
polynomial  must  be  truncated  to  a  finite  number  of  terms. 
Thus,  the  number  of  elements  in  a  structure  must  be  large 
enough  so  that  the  displacement  function  for  each  element 
closely  approximates  the  exact  displacements  in  that 
particular  region. 

In  any  numerical  method,  the  solution  should  converge 
to  the  exact  solution  as  the  size  of  the  elements  become 
small.  For  the  displacement  formulation,  it  has  been  shown 
that  under  certain  conditions  the  solution  provides  a  lower 
bound  to  the  exact  displacements  (25).  To  assure  this 
convergence  certain  conditions  must  be  satisfied.  First, 
the  displacement  function  must  be  chosen  so  that  rigid  body 
displacements  do  not  cause  straining  of  the  element. 
Second,  the  function  must  also  be  chosen  so  that  a  constant 
state  of  strain  is  obtained  as  the  element  size  approaches 
zero.  The  simplest  polynomial  function  which  satisfies 
these  two  requirements  and  also  maintains  displacement 
continuity  between  adjacent  elements  is  the  linear 
displacement  function. 

DISPLACEMENT  FUNCTION 
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Figure  A-l  shows  a  typical  triangular  element,  m,  with 
nodes  i,j,k  numbered  in  a  counter  clockwise  direction.  The 
linear  displacement  function  which  defines  the  displacements 
within  the  element  is  given  by 


(A-l) 


where  the  constants  «<£  are  determined  from  the  six  nodal 
displacements  and  nodal  coordinates  as 


bkL 


and 


(A-2) 


where  is  the 

coefficients  ,  y • 
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in-plane  area  of  the  element.  The 
and  are  given  by 
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where  r  and  z  are  coordinates  of  the  nodal  points.  The 
other  coefficients  for  the  subscripts  "j”  and  "k"  are 
obtained  by  cyclic  permutation  of  the  subscripts  i,  j,  and 
k  • 
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ELEMENT  STRAIN  FOR  AN  A XI -SYMMETRIC  GEOMETRY 

The  total  strains  at  any  point  within  an  element  are 
defined  in  terms  of  the  displacement  derivatives  as 

*/iz 

u/r 

From  eqs  A-l  to  A-6,  the  total  strains  are  written  in  terms 


£*]w 


(A— 7 ) 


of  the  nodal  displacements  and  coordinates  as 


H  =  [B]  |u] 


(A-8) 


where  Is  the  generalized  nodal  displacement  and 

[B]  = 
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where 


-  a'/r  +  hi  *■  Cc  2/r 


ELASTIC  ANALYSIS 


For  linear  elastic  and  Isotropic  materials,  the 
relationship  between  stresses  H  ,  strains  H  and  any 
Initial  strains  Is  given  by 


(A-10) 


where  t  01  is  the  elastic  material  property  matrix.  The 
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ma tr ix  [  DJ  for  axisymme tr ic  conditions  is  given  by 


S ( » -v) 

O  *v>ci-2v) 


METHOD  of  SOLUTION 


I  v/-v  O 

.zti-v),. 


(A-ll) 


The  equation  which  governs  the  elastic  response  of  a 
discretized  structure  can  be  derived  from  the  principle  of 
virtual  work  (25)  and  is  given  by 


M  H  =  |pj  4. 


where  [  K  1  is  the  elastic  stiffness  matrix  of  the  structure 
$0^13  the  generalized  displacement  vector,  n  is  the 
external  applied  load  vector,  and  1 8  the  force  vector 
due  to  the  presence  of  initial  strain. 

The  coefficients  of  the  elastic  stiffness  matrix  are 
obtained  from 


|^|  A 

=  f_Z7T\  CbJTLPJ  Cb3  rJtrJ 

M»l  J 


(A-13) 


where  the  integration  is  taken  over  the  volume  of  each 
element  and  the  summation  is  over  all  elements  in  the 
structure.  The  nodal  forces  due  to  initial  strains  are 
given  by 

ro  -  1  f  cbjtlp3  Jv.t  ( A- 1 4 ) 

^  Hi  si  J 
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Since  the  integrals  contain  terms  which  are  dependant 
on  the  radius  r,  an  efficient  method  of  integration  must  be 
chosen*  The  simplest  approximation  procedure  is  to  evaluate 
t0i  for  a  centroidal  point 

F  *  (n  *rs  *0/3  (a- 1 5 ) 

e  *  (*;  *  <A-‘6) 

In  this  case,  if  [0]  is  defined  to  be  [si  evaluated  at  the 
centroidal  point 

Ck]  a  £  z*J~Cb]t0>1Cb]  f 


f*??  ~  'L.2*^**  L&1  r 

J  Hi=  / 


(A-l 8 ) 


In  the  elastic  problem,  the  unknowns  are  the 
displacements  H  •  The  solution  for  these  unknowns  is  found 
by  the  inversion  of  Eqn  A-12* 


VISCOELASTIC  PROBLEMS 


As  discussed  in  Chapter  3,  the  displacements,  strains, 


and  stresses  are  seen  to  be  time  dependant  in  viscoelastic 
problems.  It  is  also  seen  that  the  constitutive  relation 
for  a  three-parameter  solid  can  be  integrated  in  time  to 
give  Eqn  2-50  which  repeated  here  is 


(A-19) 
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(A-20) 


If 


where  [A]  is  given  by 

\  -vU  o 

o 

o 

I  -V. 

2. 

is  treated  as  an  incremental  initial  strain 
,  then  it  is  readily  seen  that  the  viscoelastic 
constitutive  relation  enters  into  the  governing  equation 
(A-12)  only  In  the  residual  force  term  •  Then,  as  time 
is  Incremented  by  imposing  a  small  At,  the  residual  force 
term  Is  seen  to  vary  with  time.  This  means  that  Eqn  A-12 
must  be  solved  at  each  time  step  for  the  values  of 


iui ,  *$**1  >  i^j 


(A-2 1 ) 


These  values  of  M  and  f*  i  are  then  substituted  in  to 
Eqn  A-19  ,  time  is  Incremented,  and  the  process  continues 


until  a  maximum  time  is  reached.  Figure  A-2  presents  a  flow 
chart  of  the  procedure. 


APPENDIX  B 


The  application  of  the  general  correspondence  principle 
to  obtain  the  viscoelastic  solution  to  the  thick-walled 
cylinder  problem  proceeds  as  follows: 


u.  =  -e 


O 


£. 

a 


'  r  J 


(  B  —  1 ) 


Applying  the  assumption  that  the  volume  change  is  purely 
elastic,  while  the  distortion  is  a  three-parameter  solid, 
leads  to  the  following  operators: 


1  ;  -2  =  3K  ;  Q-1  (B-2) 


Substituting  into  Eqn  B-i  yields 
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P-  i^-1  3r-±lpil*t  +  < b-6 > 

y  +  y** J 
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r  libk+p)* (USffqf’  (bZ+yj+faUf+tt)*  (B-7 ) 
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At  this  point,  the  method  of  load  application  must  be 
considered*  If  the  load  is  applied  instantaneously,  the 
value  of  p  becomes  pmp/j^  *  Substituting  and  rearranging 
leads  to 


The  equation  in  this  form  is  conveniently  inverted  to  the 
time  domain  by  the  following  expressions: 


4  _ 

A  -bS- 


(B-9) 


and 
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(  B  —  1 1  ) 


u  e 

Uy  j  (,Vf+%. 
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V  L  %r  J 
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6Kp,+*\  '  '  r' 

To  obtain  the  value  of  (T  ,  consider  the  transformed 


equation 


-  jusr 
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(B-l 2 ) 


Substituting  in  the  polynomial  operators  from  Eqn  B-2  yields 


6K0  */.*) 
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( B-l 3 ) 
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1  (C>IC  +  ^Kpi+^iKJ 

Again  assuming  instantaneous  loading  leads  to 

Using  the  same  inversion  expressions  gives 

01  =  -  e*]  +.  e“*V 


(B-l 4 ) 
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where,  again 
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